# =============================================================================
# Statistical Analysis Code for:
# "Chronic adaptive versus conventional DBS response patterns in Parkinson's
#  disease: A pilot randomized crossover trial"
#
# Visualization: Sensitivity analysis results comparison plots
# Author: Jun Tanimura, MD, MSc.
# Created: 2025-11-27
# Last Updated: 2025-11-29
# Version: v2.0
#
# This script visualizes sensitivity analysis results comparing:
# 1. Main ANCOVA (full crossover model)
# 2. Reduced ANCOVA (without Period/Sequence)
# 3. Wilcoxon signed-rank test (most naive analysis)
# =============================================================================

# Load required packages
library(tidyverse)
library(ggplot2)

# Avoid conflicts by explicitly using scales functions
require(scales)

# Read the sensitivity analysis comparison results
comparison_df <- read_csv("results_sensitivity_analysis_comparison.csv")

# Define clinical evaluation categories and clean names (same as main plot)
evaluation_categories <- tribble(
  ~Evaluation, ~Category, ~Clean_Name, ~Order,
  "Mean_corr_ON_time", "Motor Diary", "ON Time", 1,
  "Mean_corr_OFF_time", "Motor Diary", "OFF Time", 2,
  "Mean_corr_Dys_sev_time", "Motor Diary", "Troublesome Dyskinesia Time", 3,
  "Mean_corr_Dys_weak_time", "Motor Diary", "Non-troublesome Dyskinesia Time", 4,
  "ADL", "Functional", "ADL", 5,
  "UPDRS_part_1", "UPDRS", "UPDRS Part I", 6,
  "UPDRS_part_2", "UPDRS", "UPDRS Part II", 7,
  "UPDRS_part_3", "UPDRS", "UPDRS Part III", 8,
  "UPDRS_part_4", "UPDRS", "UPDRS Part IV", 9,
  "UPDRS_total", "UPDRS", "UPDRS Total Score", 10,
  "MMSE", "Cognitive", "MMSE", 11,
  "PDQ-39", "QOL", "PDQ 39", 12
)

# Prepare data for plotting
plot_data <- comparison_df %>%
  left_join(evaluation_categories, by = "Evaluation") %>%
  filter(!is.na(Clean_Name)) %>%
  mutate(
    Clean_Name = factor(Clean_Name, levels = evaluation_categories$Clean_Name[evaluation_categories$Order]),
    # Direction (based on Main analysis)
    Direction = ifelse(Main_Cohens_D > 0, "aDBS_Superior", "cDBS_Superior"),
    Direction = factor(Direction, levels = c("aDBS_Superior", "cDBS_Superior"))
  ) %>%
  arrange(Order)

# Create long format for plotting
plot_data_combined <- bind_rows(
  plot_data %>%
    select(Evaluation, Clean_Name, Order, Direction, 
           Effect_Size = Main_Cohens_D, 
           CI_Lower = Main_Cohens_D_CI_Lower,
           CI_Upper = Main_Cohens_D_CI_Upper,
           P = Main_P) %>%
    mutate(Analysis_Type = "Main\n(Full Crossover)"),
  
  plot_data %>%
    select(Evaluation, Clean_Name, Order, Direction,
           Effect_Size = Reduced_Cohens_D,
           CI_Lower = Reduced_Cohens_D_CI_Lower,
           CI_Upper = Reduced_Cohens_D_CI_Upper,
           P = Reduced_P) %>%
    mutate(Analysis_Type = "Reduced\n(No Period/Seq)"),
  
  plot_data %>%
    select(Evaluation, Clean_Name, Order, Direction,
           Effect_Size = Wilcoxon_Cohens_D,
           CI_Lower = Wilcoxon_Cohens_D_CI_Lower,
           CI_Upper = Wilcoxon_Cohens_D_CI_Upper,
           P = Wilcoxon_P) %>%
    mutate(Analysis_Type = "Wilcoxon\n(Naive)")
) %>%
  mutate(
    Analysis_Type = factor(Analysis_Type, 
                          levels = c("Main\n(Full Crossover)", 
                                     "Reduced\n(No Period/Seq)", 
                                     "Wilcoxon\n(Naive)")),
    # Order for plotting (bottom to top)
    Clean_Name = factor(Clean_Name, levels = rev(evaluation_categories$Clean_Name[evaluation_categories$Order]))
  )

# Analysis type colors (purple, gray, yellow scheme - muted tones)
analysis_colors <- c(
  "Main\n(Full Crossover)" = "#8B7AB8",      # Muted purple
  "Reduced\n(No Period/Seq)" = "#808080",    # Medium gray
  "Wilcoxon\n(Naive)" = "#C9A961"            # Muted golden yellow
)

# ===========================================
# MAIN FIGURE: Comparison Forest Plot
# ===========================================

# Forest plot for sensitivity analysis comparison on effect sizes (Cohen's d) with 95% CI.
create_sensitivity_comparison_plot <- function(data) {
  
  # Calculate dynamic x-axis limits
  min_value <- min(c(data$Effect_Size, data$CI_Lower), na.rm = TRUE)
  max_value <- max(c(data$Effect_Size, data$CI_Upper), na.rm = TRUE)
  x_range <- max_value - min_value
  x_limits <- c(min_value - x_range * 0.15, max_value + x_range * 0.15)
  
  # Create the comparison plot (y-axis remains as factor, use position_dodge for multiple points)
  p <- ggplot(data, aes(y = Clean_Name, x = Effect_Size, color = Analysis_Type)) +
    
    # Add vertical line at zero
    geom_vline(xintercept = 0, linetype = "solid", color = "gray50", linewidth = 0.5) +
    
    # Add dotted lines for effect size thresholds
    {if(x_limits[2] >= 0.2) geom_vline(xintercept = c(-0.2, 0.2), linetype = "dotted", color = "gray70", alpha = 0.7)} +
    {if(x_limits[2] >= 0.5) geom_vline(xintercept = c(-0.5, 0.5), linetype = "dotted", color = "gray70", alpha = 0.7)} +
    {if(x_limits[2] >= 0.8) geom_vline(xintercept = c(-0.8, 0.8), linetype = "dotted", color = "gray70", alpha = 0.7)} +
    
    # Add error bars for 95% CI (simple horizontal lines without caps)
    # Using geom_errorbarh with height = 0 to remove caps, position_dodge for spacing
    geom_errorbarh(aes(xmin = CI_Lower, xmax = CI_Upper),
                   height = 0, linewidth = 0.9, alpha = 0.5,
                   position = position_dodge(width = 0.5)) +
    
    # Add points for effect sizes (all circles, differentiated by color gradient)
    geom_point(shape = 19,  # filled circle for all
               size = 3, alpha = 0.8,
               position = position_dodge(width = 0.5)) +
    
    # Color scale (gray gradient)
    scale_color_manual(
      values = analysis_colors,
      name = "Analysis Type",
      guide = guide_legend(title.position = "top", title.hjust = 0.5)
    ) +
    
    # Ensure y-axis is discrete (factor labels, not numeric)
    scale_y_discrete(drop = FALSE) +
    
    # Dynamic X-axis
    scale_x_continuous(
      limits = x_limits,
      breaks = scales::pretty_breaks(n = 6),
      labels = function(x) sprintf("%.1f", x)
    ) +
    
    # Labels
    labs(
      title = "Sensitivity Analysis: Effect Sizes Comparison",
      subtitle = "Main vs Reduced ANCOVA vs Wilcoxon Signed-Rank",
      x = "Cohen's d (95% CI)",
      y = "",
      caption = "Positive values indicate aDBS superiority\nError bars represent 95% confidence intervals\nMain: Full crossover ANCOVA | Reduced: No Period/Sequence | Wilcoxon: Naive paired test"
    ) +
    
    # Theme
    theme_classic(base_size = 12) +
    theme(
      plot.title = element_text(size = 14, face = "bold"),
      plot.subtitle = element_text(size = 11, color = "gray40"),
      axis.text.y = element_text(size = 9),
      axis.text.x = element_text(size = 9),
      axis.title.x = element_text(size = 10),
      panel.grid.major.x = element_blank(),
      plot.caption = element_text(size = 8, color = "gray50", hjust = 0),
      plot.margin = margin(10, 10, 10, 10),
      legend.position = "bottom",
      legend.text = element_text(size = 9),
      legend.title = element_text(size = 10, face = "bold")
    )
  
  return(p)
}

# ===========================================
# CREATE AND SAVE PLOTS
# ===========================================

# Generate main comparison plot
main_comparison_plot <- create_sensitivity_comparison_plot(plot_data_combined)

# Save plot (narrower and taller for portrait orientation)
ggsave("figure_sensitivity_comparison_forest_plot.pdf", 
       main_comparison_plot, width = 4, height = 5.5, bg = "white")

# Print summary
cat("\n=== SENSITIVITY ANALYSIS VISUALIZATION COMPLETE ===\n")
cat("Files created:\n")
cat("1. figure_sensitivity_comparison_forest_plot.pdf - Combined comparison plot\n")

# Display main plot
print(main_comparison_plot)

cat("\nVisualization completed successfully!\n")

# Clean up (scales package cannot be unloaded as it's imported by ggplot2)
# No cleanup needed

